#include<iostream>
#include"rpa.h"
#include"constants.h"
#include"help.h"
RPA::RPA(const System_Tz_HF* Ptr):pSystem_Tz_HF(Ptr),isTDA(false)
{
  vector<TwoBodyState_ph> temp;
  Jmax=0;
  while(1)
    {
      if(set_ph_same(Jmax,0,temp)>0||set_ph_same(Jmax,1,temp)>0
	 ||set_ph_np(Jmax,0,temp)>0||set_ph_np(Jmax,1,temp)>0
	 ||set_ph_pn(Jmax,0,temp)>0||set_ph_pn(Jmax,1,temp)>0
	 )
	++Jmax;
      else
	break;
    }
  ph_sames.resize(2*Jmax);
  ph_nps.resize(2*Jmax);
  ph_pns.resize(2*Jmax);

  A_sames.resize(2*Jmax);
  B_sames.resize(2*Jmax);
  X_sames.resize(2*Jmax);
  Y_sames.resize(2*Jmax);
  Omega_sames.resize(2*Jmax);

  Anpnps.resize(2*Jmax);
  Apnpns.resize(2*Jmax);
  Bnppns.resize(2*Jmax);
  T_plus_Xs.resize(2*Jmax);
  T_plus_Ys.resize(2*Jmax);
  T_plus_Omegas.resize(2*Jmax);
  T_minus_Xs.resize(2*Jmax);
  T_minus_Ys.resize(2*Jmax);
  T_minus_Omegas.resize(2*Jmax);
}
int RPA::set_ph_same(int J,int Par,vector<TwoBodyState_ph>&ph_same)
{
  const  HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_same.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
      	  if(HFOrbs[h].tz!=HFOrbs[p].tz) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_same.push_back(temp);
	}
    }
  return ph_same.size();
}
int RPA::set_ph_np(int J,int Par,vector<TwoBodyState_ph>&ph_np)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_np.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
	  //tz=-1 for proton,1 for neutron
	  if((HFOrbs[h].tz!=-1)||(HFOrbs[p].tz!=1)) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_np.push_back(temp);
	}
    }
  return ph_np.size();
}
int RPA::set_ph_pn(int J,int Par,vector<TwoBodyState_ph>&ph_pn)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  ph_pn.clear();
  for(int h=0;h<=pSystem_Tz_HF->FermiSurface;++h)
    {
      for(int p=pSystem_Tz_HF->FermiSurface+1;p<HFOrbs.size();++p)
	{
	  //tz=-1 for proton,1 for neutron
	  if((HFOrbs[h].tz!=1)||(HFOrbs[p].tz!=-1)) continue;
	  if((HFOrbs[h].l+HFOrbs[p].l-Par)%2!=0) continue;
	  if(abs(HFOrbs[h].jj-HFOrbs[p].jj)>2*J) continue;
	  if((HFOrbs[h].jj+HFOrbs[p].jj)<2*J) continue;
	  TwoBodyState_ph temp(p,h);
	  ph_pn.push_back(temp);
	}
    }
  return ph_pn.size();
}

double RPA::a(int bra_p,int bra_h,int ket_p,int ket_h,int J)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  double temp=0.0;
  if((bra_h==ket_h)&&(bra_p==ket_p))
    temp=HFOrbs[bra_p].e-HFOrbs[bra_h].e;
  int bra_p_jj=HFOrbs[bra_p].jj;int bra_h_jj=HFOrbs[bra_h].jj;
  int ket_p_jj=HFOrbs[ket_p].jj;int ket_h_jj=HFOrbs[ket_h].jj;
  int j_min=max(abs(ket_p_jj-bra_h_jj),abs(bra_p_jj-ket_h_jj))/2;
  int j_max=min(ket_p_jj+bra_h_jj,bra_p_jj+ket_h_jj)/2;
  for(int k=j_min;k<=j_max;++k)
    {
      int phase=((ket_p_jj+bra_h_jj)/2-k)%2?-1:1;
      temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_p_jj,ket_h_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,ket_h,bra_h,ket_p,k);
    }
  return temp;
}

double RPA::b(int bra_p,int bra_h,int ket_p,int ket_h,int J)
{
  const HFOrbitals & HFOrbs=pSystem_Tz_HF->Orbitals;
  double temp=0.0;
  int bra_p_jj=HFOrbs[bra_p].jj;int bra_h_jj=HFOrbs[bra_h].jj;
  int ket_p_jj=HFOrbs[ket_p].jj;int ket_h_jj=HFOrbs[ket_h].jj;
  int j_min=max(abs(ket_p_jj-bra_p_jj),abs(bra_h_jj-ket_h_jj))/2;
  int j_max=min(ket_p_jj+bra_p_jj,bra_h_jj+ket_h_jj)/2;
  for(int k=j_min;k<=j_max;++k)
    {
      int phase=((ket_p_jj+bra_h_jj)/2+k+J)%2?-1:1;
      temp+=phase*(2*k+1)*gsl_sf_coupling_6j(bra_p_jj,bra_h_jj,2*J,ket_h_jj,ket_p_jj,2*k)*pSystem_Tz_HF->get2B(bra_p,ket_p,bra_h,ket_h,k);
    }
  return temp;
}

void RPA::set_AB_same(int J,int Par)
{
  set_ph_same(J,Par,ph_sames[2*J+Par]);
  int dim=ph_sames[2*J+Par].size();
  A_sames[2*J+Par].resize(dim,dim);
  B_sames[2*J+Par].resize(dim,dim);
  B_sames[2*J+Par]=MatrixXd::Zero(dim,dim);

  for(int i=0;i<dim;++i)
    {
      for(int j=i;j<dim;++j)
	{
	  int bra_p=ph_sames[2*J+Par][i].a;
	  int bra_h=ph_sames[2*J+Par][i].b;
	  int ket_p=ph_sames[2*J+Par][j].a;
	  int ket_h=ph_sames[2*J+Par][j].b;
	  A_sames[2*J+Par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
	  A_sames[2*J+Par](j,i)=A_sames[2*J+Par](i,j);
	  if(!isTDA)
	    {
	      B_sames[2*J+Par](i,j)=b(bra_p,bra_h,ket_p,ket_h,J);
	      B_sames[2*J+Par](j,i)=B_sames[2*J+Par](i,j);
	    }
	}
    }
}

void RPA::cal_same(int J,int Par)
{
  set_AB_same(J,Par);
  if(ph_sames[2*J+Par].empty())
    {
      A_sames[2*J+Par].resize(0,0);
      B_sames[2*J+Par].resize(0,0);
      X_sames[2*J+Par].resize(0,1);
      Y_sames[2*J+Par].resize(0,1);
      Omega_sames[2*J+Par].resize(0,1);
      return;
    }
  
  int dim=ph_sames[2*J+Par].size();
  MatrixXcd R(2*dim,2*dim);
  for(int i=0;i<dim;++i)
    {
      for(int j=0;j<dim;++j)
	{
	  R(i,j)=complexd(A_sames[2*J+Par](i,j),0);
	  R(dim+i,dim+j)=complexd(-A_sames[2*J+Par](i,j),0);
	  R(i,dim+j)=complexd(B_sames[2*J+Par](i,j),0);
	  R(dim+i,j)=complexd(-B_sames[2*J+Par](i,j),0);
	}
    }
  ComplexEigenSolver<MatrixXcd> es(R);
  VectorXcd eigen_values=es.eigenvalues();
  MatrixXcd eigen_vectors=es.eigenvectors();
  Omega_sames[2*J+Par].resize(dim,1);
  X_sames[2*J+Par].resize(dim,dim);
  Y_sames[2*J+Par].resize(dim,dim);
  int j=0;
  //  int zero_count=0;
  for(int i=0;i<2*dim;++i)
    {
      double norm=eigen_vectors.col(i).head(dim).real().squaredNorm()-eigen_vectors.col(i).tail(dim).real().squaredNorm();
      if(norm>0)
      	{
      	  X_sames[2*J+Par].col(j)=eigen_vectors.col(i).head(dim).real();
      	  Y_sames[2*J+Par].col(j)=eigen_vectors.col(i).tail(dim).real();
      	  X_sames[2*J+Par].col(j)/=sqrt(norm);
      	  Y_sames[2*J+Par].col(j)/=sqrt(norm);
      	  Omega_sames[2*J+Par](j)=eigen_values(i).real();
      	  j++;
      	}
      // //to handle the spurious state of negative 1 caused by the breaking symmetry of translation.
      // if(abs(eigen_values(i).real())<1)
      // 	{
      // 	  if(zero_count)
      // 	    continue;
      // 	  else
      // 	    {
      // 	      ++zero_count;
      // 	      Omega_sames[2*J+Par](j)=eigen_values(i).real();
      // 	      X_sames[2*J+Par].col(j)=eigen_vectors.col(i).head(dim).real();
      // 	      Y_sames[2*J+Par].col(j)=eigen_vectors.col(i).tail(dim).real();
      // 	      double norm=X_sames[2*J+Par].col(j).squaredNorm()-Y_sames[2*J+Par].col(j).squaredNorm();
      // 	      if(norm>0)
      // 		{
      // 		  X_sames[2*J+Par].col(j)/=sqrt(norm);
      // 		  Y_sames[2*J+Par].col(j)/=sqrt(norm);
      // 		}
      // 	      else
      // 		{
      // 		  Omega_sames[2*J+Par](j)=0;
      // 		}
      // 	      ++j;

      // 	    }
      // 	}
      // else if(eigen_values(i).real()>0)
      // 	{
      // 	  Omega_sames[2*J+Par](j)=eigen_values(i).real();
      // 	  X_sames[2*J+Par].col(j)=eigen_vectors.col(i).head(dim).real();
      // 	  Y_sames[2*J+Par].col(j)=eigen_vectors.col(i).tail(dim).real();
      // 	  double norm=X_sames[2*J+Par].col(j).squaredNorm()-Y_sames[2*J+Par].col(j).squaredNorm();
      // 	  if(norm>0)
      // 	    {
      // 	      X_sames[2*J+Par].col(j)/=sqrt(norm);
      // 	      Y_sames[2*J+Par].col(j)/=sqrt(norm);
      // 	    }
      // 	  else
      // 	    {
      // 	      Omega_sames[2*J+Par](j)=0;
      // 	    }
      // 	  ++j;
      // 	}
    }
}


void RPA::set_AB(int J,int Par)
{
  set_ph_np(J,Par,ph_nps[2*J+Par]);
  set_ph_pn(J,Par,ph_pns[2*J+Par]);
  int dim_np=ph_nps[2*J+Par].size();
  int dim_pn=ph_pns[2*J+Par].size();
   
  Anpnps[2*J+Par].resize(dim_np,dim_np);
  Apnpns[2*J+Par].resize(dim_pn,dim_pn);
  Bnppns[2*J+Par].resize(dim_np,dim_pn);

  //setup Anpnp
  for(int i=0;i<dim_np;++i)
    {
      for(int j=i;j<dim_np;++j)
	{
	  int bra_p=ph_nps[2*J+Par][i].a;
	  int bra_h=ph_nps[2*J+Par][i].b;
	  int ket_p=ph_nps[2*J+Par][j].a;
	  int ket_h=ph_nps[2*J+Par][j].b;
	  Anpnps[2*J+Par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
	  Anpnps[2*J+Par](j,i)=Anpnps[2*J+Par](i,j);
	}
    }

  
  //setup Apnpn
  for(int i=0;i<dim_pn;++i)
    {
      for(int j=i;j<dim_pn;++j)
	{
	  int bra_p=ph_pns[2*J+Par][i].a;
	  int bra_h=ph_pns[2*J+Par][i].b;
	  int ket_p=ph_pns[2*J+Par][j].a;
	  int ket_h=ph_pns[2*J+Par][j].b;
	  Apnpns[2*J+Par](i,j)=a(bra_p,bra_h,ket_p,ket_h,J);
	  Apnpns[2*J+Par](j,i)=Apnpns[2*J+Par](i,j);
	}
    }

  if(isTDA)
    {
      Bnppns[2*J+Par]=MatrixXd::Zero(dim_np,dim_pn);
      return;
    }

  //setup Bnppn
  for(int i=0;i<dim_np;++i)
    {
      for(int j=0;j<dim_pn;++j)
	{
	  int bra_p=ph_nps[2*J+Par][i].a;
	  int bra_h=ph_nps[2*J+Par][i].b;
	  int ket_p=ph_pns[2*J+Par][j].a;
	  int ket_h=ph_pns[2*J+Par][j].b;
	  Bnppns[2*J+Par](i,j)=b(bra_p,bra_h,ket_p,ket_h,J);
	}
    }
}

void RPA::cal_ce(int J,int Par)
{
  set_AB(J,Par);
  if(ph_nps[2*J+Par].empty()&&ph_pns[2*J+Par].empty())
    {
      Anpnps[2*J+Par].resize(0,0);
      Apnpns[2*J+Par].resize(0,0);
      Bnppns[2*J+Par].resize(0,0);
      T_plus_Xs[2*J+Par].resize(0,1);
      T_plus_Ys[2*J+Par].resize(0,1);
      T_plus_Omegas[2*J+Par].resize(0,1);

      T_minus_Xs[2*J+Par].resize(0,1);
      T_minus_Ys[2*J+Par].resize(0,1);
      T_minus_Omegas[2*J+Par].resize(0,1);

      return;
    }
  int dim_np=ph_nps[2*J+Par].size();
  int dim_pn=ph_pns[2*J+Par].size();

  MatrixXcd R(dim_np+dim_pn,dim_np+dim_pn);
  for(int i=0;i<dim_np;++i)
    {
      for(int j=0;j<dim_np;++j)
  	{
  	  R(i,j)=complexd(Anpnps[2*J+Par](i,j),0);
  	}
      for(int j=0;j<dim_pn;++j)
	{
	  R(i,dim_np+j)=complexd(Bnppns[2*J+Par](i,j),0);
	}
    }
  for(int i=0;i<dim_pn;++i)
    {
      for(int j=0;j<dim_np;++j)
  	{
  	  R(dim_np+i,j)=complexd(-Bnppns[2*J+Par](j,i),0);
  	}
      for(int j=0;j<dim_pn;++j)
	{
	  R(dim_np+i,dim_np+j)=complexd(-Apnpns[2*J+Par](i,j),0);
	}
    }
  ComplexEigenSolver<MatrixXcd> es(R);
  VectorXcd eigen_values=es.eigenvalues();
  MatrixXcd eigen_vectors=es.eigenvectors();
  int T_plus_dim=0,T_minus_dim=0;
  for(int i=0;i<eigen_values.size();i++)
    {
      double temp1=eigen_vectors.col(i).head(dim_np).real().squaredNorm();
      double temp2=eigen_vectors.col(i).tail(dim_pn).real().squaredNorm();
      if(temp1>temp2)
	++T_plus_dim;
      else
	++T_minus_dim;
    }
  T_plus_Omegas[2*J+Par].resize(T_plus_dim,1);
  T_minus_Omegas[2*J+Par].resize(T_minus_dim,1);

  T_plus_Xs[2*J+Par].resize(dim_np,T_plus_dim);
  T_plus_Ys[2*J+Par].resize(dim_pn,T_plus_dim);

  T_minus_Ys[2*J+Par].resize(dim_np,T_minus_dim);
  T_minus_Xs[2*J+Par].resize(dim_pn,T_minus_dim);

  int i_plus=0,i_minus=0;
  for(int i=0;i<eigen_values.size();++i)
    {
      double temp1=eigen_vectors.col(i).head(dim_np).real().squaredNorm();
      double temp2=eigen_vectors.col(i).tail(dim_pn).real().squaredNorm();
      double norm=0;
      if(temp1>temp2)
	{
	  norm=sqrt(temp1-temp2);
	  T_plus_Omegas[2*J+Par](i_plus)=eigen_values(i).real();
	  T_plus_Xs[2*J+Par].col(i_plus)=eigen_vectors.col(i).head(dim_np).real();
	  T_plus_Ys[2*J+Par].col(i_plus)=eigen_vectors.col(i).tail(dim_pn).real();
	  T_plus_Xs[2*J+Par].col(i_plus)/=norm;
	  T_plus_Ys[2*J+Par].col(i_plus)/=norm;
	  ++i_plus;
	}
      else
	{
	  norm=sqrt(temp2-temp1);
	  T_minus_Omegas[2*J+Par](i_minus)=-eigen_values(i).real();
	  T_minus_Ys[2*J+Par].col(i_minus)=eigen_vectors.col(i).head(dim_np).real();
	  T_minus_Xs[2*J+Par].col(i_minus)=eigen_vectors.col(i).tail(dim_pn).real();
	  T_minus_Ys[2*J+Par].col(i_minus)/=norm;
	  T_minus_Xs[2*J+Par].col(i_minus)/=norm;
	  ++i_minus;
	}
    }
  return;
}
void RPA::cal(int type,int J,int Par)
{
  type?cal_ce(J,Par):cal_same(J,Par);
}
void RPA::cal()
{
  for(int J=0;J<Jmax;J++)
    {
      for(int Par=0;Par<2;Par++)
	{
	  for(int Type=0;Type<2;Type++)
	    {
	      cal(Type,J,Par);
	    }
	}
    }
}

void RPA::setTDAflag()
{
  isTDA=true;
}

double RPA::getDeltaE()
{
  //two ways to cal. correlation energy as a check
  double e_rpa1=0,e_rpa2=0;
  double Etotal1=0,Etotal2=0;
  for(int J=0;J<Jmax;++J)
    {
      e_rpa1=0;
      e_rpa2=0;
      for(int par=0;par<2;par++)
	{
	  for(int i=0;i<Omega_sames[2*J+par].size();++i)
	    {
	      e_rpa1-=Omega_sames[2*J+par](i)*Y_sames[2*J+par].col(i).squaredNorm();
	      e_rpa2+=0.5*( Omega_sames[2*J+par](i)-A_sames[2*J+par](i,i) );
	    }
	}
      Etotal1+=(2*J+1)*e_rpa1;
      Etotal2+=(2*J+1)*e_rpa2;
      //cout<<J<<"\t"<<e_rpa1<<"\t"<<e_rpa2<<endl;
    }
  for(int J=0;J<Jmax;++J)
    {
      e_rpa1=0;
      e_rpa2=0;
      for(int par=0;par<2;par++)
	{
	  for(int i=0;i<T_plus_Omegas[2*J+par].size();++i)
	    {
	      e_rpa1-=T_plus_Omegas[2*J+par](i)*T_plus_Ys[2*J+par].col(i).squaredNorm();
	      e_rpa2+=0.5*T_plus_Omegas[2*J+par](i);
	    }
	  for(int i=0;i<ph_nps[2*J+par].size();i++)
	    {
	      e_rpa2+=-0.5*Anpnps[2*J+par](i,i);
	    }
	  for(int i=0;i<T_minus_Omegas[2*J+par].size();++i)
	    {
	      e_rpa1+=-T_minus_Omegas[2*J+par](i)*T_minus_Ys[2*J+par].col(i).squaredNorm();
	      e_rpa2+=0.5*T_minus_Omegas[2*J+par](i);
	    }
	  for(int i=0;i<ph_pns[2*J+par].size();i++)
	    {
	      e_rpa2+=-0.5*Apnpns[2*J+par](i,i);
	    }
	  //	  cout<<J<<"\t"<<e_rpa1<<"\t"<<e_rpa2<<endl;
	}
      Etotal1+=(2*J+1)*e_rpa1;
      Etotal2+=(2*J+1)*e_rpa2;
    }
  //  cout<<"Results of correlation energy using two formulas:\n";
  //cout<<Etotal1<<"\t"<<Etotal2<<endl;
  return Etotal1;
}
double RPA::calOccuNum()
{
  OccuNum.resize(pSystem_Tz_HF->Orbitals.size());
  for(int i=0;i<pSystem_Tz_HF->Orbitals.size();i++)
    {
      OccuNum[i]=0;
    }
  for(int J=0;J<Jmax;J++)
    {
      for(int Par=0;Par<2;Par++)
	{
	  for(int i=0;i<Y_sames[2*J+Par].cols();++i)
	    {
	      for(int j=0;j<ph_sames[2*J+Par].size();j++)
	  	{
	  	  int p=ph_sames[2*J+Par][j].a;
	  	  int h=ph_sames[2*J+Par][j].b;
	  	  OccuNum[p]+=(2*J+1)*pow(Y_sames[2*J+Par](j,i),2);
	  	  OccuNum[h]+=(2*J+1)*pow(Y_sames[2*J+Par](j,i),2);
	  	}
	    }
	  for(int i=0;i<T_plus_Ys[2*J+Par].cols();i++)
	    {
	      for(int j=0;j<ph_pns[2*J+Par].size();j++)
	  	{
	  	  int p=ph_pns[2*J+Par][j].a;
	  	  int h=ph_pns[2*J+Par][j].b;
	  	  OccuNum[p]+=(2*J+1)*pow(T_plus_Ys[2*J+Par](j,i),2);
	  	  OccuNum[h]+=(2*J+1)*pow(T_plus_Ys[2*J+Par](j,i),2);
	  	}
	    }
	  for(int i=0;i<T_minus_Ys[2*J+Par].cols();i++)
	    {
	      for(int j=0;j<ph_nps[2*J+Par].size();j++)
	  	{
	  	  int p=ph_nps[2*J+Par][j].a;
	  	  int h=ph_nps[2*J+Par][j].b;
	  	  OccuNum[p]+=(2*J+1)*pow(T_minus_Ys[2*J+Par](j,i),2);
	  	  OccuNum[h]+=(2*J+1)*pow(T_minus_Ys[2*J+Par](j,i),2);
	  	}
	    }
	}
    }
  for(int orb=0;orb<pSystem_Tz_HF->Orbitals.size();orb++)
    {
      if(orb>pSystem_Tz_HF->FermiSurface)
	{
	  OccuNum[orb]*=0.5;
	}
      else
	{
	  OccuNum[orb]=(pSystem_Tz_HF->Orbitals[orb].jj+1)-0.5*OccuNum[orb];
	}
    }
}
void RPA::set_q(int J,int Par,int T)
{
  double Rms=pSystem_Tz_HF->rms();
  int dim=ph_sames[2*J+Par].size();
  q.resize(dim);
  if(J==1)
    { 
      if(T==0)
	{
	  for(int i=0;i<dim;++i)
	    {
	      q[i]=pSystem_Tz_HF->Q1_s(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b,Rms);
	    }
	}
      else
	{
	  for(int i=0;i<dim;++i)
	    {
	      q[i]=pSystem_Tz_HF->Q1_v(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b);
	    }
	}
    }
  else
    {
      for(int i=0;i<dim;++i)
	{
	  q[i]=pSystem_Tz_HF->Q(ph_sames[2*J+Par][i].a,ph_sames[2*J+Par][i].b,J,T);
	}
    }
}
double RPA::BEJ(int J,int Par,int v)
{
  double temp=0;
  int dim=ph_sames[2*J+Par].size();
  for(int i=0;i<dim;++i)
    {
      temp+=(X_sames[2*J+Par](i,v)+phase(J)*Y_sames[2*J+Par](i,v))*q[i];
    }
  temp=temp*temp;
  return temp;
}
double RPA::S(int J,int Par,int T)
{
  set_q(J,Par,T);
  double temp=0;
  for(int i=0;i<Omega_sames[2*J+Par].size();++i)
    {
      temp+=Omega_sames[2*J+Par](i)*BEJ(J,Par,i);
    }
  return temp;
}
void RPA::R(int J,int Par,int T,ostream&fout,const double gamma,int points,const double dE)
{
  set_q(J,Par,T);
  double E=0;
  for(int i=0;i<points;++i)
    {
      E=i*dE;
      double temp=0;
      for(int j=0;j<Omega_sames[2*J+Par].size();++j)
	{
	  temp+=0.5*gamma/( (E-Omega_sames[2*J+Par](j))*(E-Omega_sames[2*J+Par](j)) + 0.25*gamma*gamma )/Pi * BEJ(J,Par,j);
	}
      fout<<E<<"\t"<<temp<<endl;
    }
}
